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ABSTRAGT 

Large scale structure measurements require accurate and precise knowledge of the survey depth 
— typically expressed in the form of a limiting magnitude — as a function of position on the sky. 
To date, most surveys only compute the point-source limiting magnitude measured within a fixed 
metric aperture. However, this quantity is ill suited to describe the limiting depth of galaxies, which 
depends on the detailed interplay of survey systematics with galaxy shapes and sizes. We describe 
an empirical method for directly estimating the limiting magnitude for large photometric surveys, 
and apply it to ~ 10,000 deg^ of SDSS DR8 data. Gombined with deeper imaging from SDSS Stripe 
82 and GFHTLens, we are able to use these depth maps to estimate the location-dependent galaxy 
detection completeness at any point within the full BOSS DR8 survey region. We show that these 
maps can be used to construct random points suitable for unbiased estimation of correlation functions 
for galaxies near the survey limiting magnitude. Finally, we provide limiting magnitude maps for 
galaxies in SDSS DR8 in HEALPix format with NSIDE=2048. 

Subject headings: surveys, galaxies: general 


1. INTRODUGTION 

In recent years, observational astronomy has been rev¬ 
olutionized by wide-held optical surveys. The impact 
of the data frorn the Sloan Digital Sky Survey (SDSS 
York et al. 2000) has been enormous, covering a range 
ot topics too broad to mention. The next era of large 
photometric surve ys is now upon us, with the D ark En¬ 
ergy Survey (DES [The DES Collaboration||2005 ) under¬ 
way, imaging 5,000 deg^ with a depth that is two mag¬ 
nitudes deeper than that of S DSS. Similar surveys like 
the Kilo-Degree Survey (KiDS |de Jong et al. 20151 and 
the Hyper-Suprime Camerrn (HSG) are also underway. 
By the end of t he decade the Larg e Synoptic Survey 
Telescope (LSST |Ivezic et al.| 12008) will start to take 
data, with unprecedented deptn over fully half the sky. 
However, as hard as we try to build homogeneous sur¬ 
veys with uniform depth, there will always be variations 
in sky coverage, with exposure time, seeing, sky bright¬ 
ness, and other factors varying over the course of nights 
and seasons. These systematic variations can impact the 
scientific products derived from these surveys, affecting 
spectroscopic target selection, optical cluster measure¬ 
ments, and correlation function measurements (to name 
but a few). 

One obvious option to mitigate these factors is to limit 
any object selection to those regions where the imag¬ 
ing is nearly complete. This of course raise the question 
of choosing which regions have high completeness. One 
option is to mask out regions which ha ve anomalously 
low object density (e.g. Liu et al. 2015). Other options 
include cross-correlating galaxy densities with known 
sources of systematic errors in order to choose which re- 
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gions were contaminated and should be masked, as for 


early studies with SDSS (Scranton et al 
Baryon Osci llation Spectroscopic Survey 


ri|200l). 
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et al. 2009) made use of relatively faint targets, re - 


quirmg more sophisticated analyses (Ross et al. 2011). 
Galaxy detection probabilities were estimated via cross¬ 
correlation with systematics maps, and galaxies used in 
the correlation function measurements were reweighted 
accordingly. However, these effects were mostly pertur¬ 
bative for the spectroscopic targets c onsidered. _ 

The WiggleZ spectroscopic survey (|Blake et al.||2010| 
had a much more challenging task due to the combination 
of SDSS and Galaxy Evolution Exploreij^ (GALEX) tar¬ 
get selection. Density variations as a function of Galac¬ 
tic dust and GALEX exposure time were modeled via a 
simple completeness function. These models were then 
used to compute the full survey selection function. Most 
recently, a significant amount of effort has gone into sim¬ 
ulating the selection function directly over a large sur¬ 
vey by inserting and measuring fake galaxies with k nown 
properties (i.e. the Balrog code,|Suchyta et al.|2015|). Al¬ 
though this may be the most robust method, it suffers 
from a couple limitations. First, it can be expensiv^ 
Second, inserting fake galaxies will tend to correlate ex¬ 
isting structure with the selection function, somewhat 
muddying the interpretation of the resulting random cat¬ 
alogs. Of course, this effect is in the real data as well, but 
since the interpretations are non-trivial, having a com¬ 
pletely independent approach for characterizing selection 
functions can be highly complementary. 

In this paper, we suggest an alternative, pragmatic ap¬ 
proach. The main idea behind our measurement is sim¬ 
ple: the magnitude error Am of a galaxy should be a 
function of the galaxy’s magnitude and the survey depth: 
the deeper the survey, the smaller the photometric error 
Urn- Given a model for magnitude errors as a function 
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of magnitude and survey depth Gm = we use 

observations of galaxies (to,Ato) over a patch of sky to 
recover the limiting magnitude TOiim of the survey patch. 
In practice, this direct measurement of the survey depth 
can only be done in relatively coarse pixels. We over¬ 
come this difficulty by fitting the recovered depth map 
as a function of the systematics map for the survey (e.g. 
sky noise, PSF size, etc), and then utilize these fine struc¬ 
ture maps to generate a final high resolution estimate of 
the survey depth. 

The layout of this paper is as follows. In Section 
we describe the SDSS data used in this paper, although 
our method is fully generalizable to any survey data. In 
Section we describe our method for measuring the sur¬ 
vey depth in coarse pixels on the sky, and in Section 
we show how we can use machine learning methods to 
reconstruct the survey depth at high resolution by mak¬ 
ing use of detailed maps of survey systematics. In Sec¬ 
tion we use deeper data to model the galaxy detection 
com^eteness, and show how this is a simple function of 
survey depth. In Section we show a practical appli¬ 
cation of our method in a simulated computation of the 
correlation function near the survey limiting magnitude 
where depth and completeness variations are signihcant. 
Finally, in Section we summarize our findings. 

2. DATA 

As discussed above, our method is designed to han¬ 
dle an arbitrary photometric galaxy catalog. As a case 
study, in this paper we determine the depth of SDSS DR8 
photometric data. 


2.1. SDSS DR8 Photometry 

The input gala xy catalog for this w ork is derived from 
SDSS DR8 data ( Aihara et al.||2011 ). This data release 
includes more than 14,000 deg^ of drift-scan imaging in 
the Northern and Southern Galactic caps. The survey 
edge us ed is the same as th at used for BOSS target se¬ 
lection ( Dawson et al.|20T3 ), which reduces the total area 
to « 10,500 deg^ with high-quality observations and a 
well-defined contiguous footprinifT 
The galaxy selection employea is nearly the same as 
that used in the construction of the SDSS DR8 redMaP- 


Per cluster catalog (Rykoff et al. 2014). The primary 
difference is that we limit the catalog to icmod < 22.0, 
rather than Zcmod < 21, to ensure that we have galax¬ 
ies that are fainter than the lOcr limiting magnitude over 
the full survey footprint. Furthermore, we have not cor¬ 
rected for Galactic reddening. We then filter all ob¬ 
jects with any of the following flags set in the g, r, or 
i bands: SATUR_CENTER, BRIGHT, TQD_MANY_PEAKS, and 
(NOT BLENDED OR NODEBLEND). 

In this work, we reconstruct the SDSS DR8 MDDEL_MAG 
limiting magnitudes in u, g, r, i and z, and the com¬ 
posite CM0DEL_MAG limiting magnitudes in r and i. The 
MODEL magnitudes are estimated by fitting each of an 
exponential and deVaucouleurs profile to each galaxy in 
the r band, and applying the better-fit model to each in¬ 
dividual band. As such, MODEL magnitudes are suitable 


for galaxy color measurements. The composite CMODEL 
magnitudes are computed as a linear combination of the 
best-fit exponential and deVaucouleurs models, and are 
useful for galaxy total magnitude measurements. In this 
paper, we denote MODEL and CMODEL magnitudes with 
subscripts “mod” and “cmod” respectively. 

2.2. SDSS Systematics Maps 

The SDSS survey consists of drift-scan imaging of a 
large number of stripe s that extend along great cir¬ 
cles (|York et al.||2000D. Because of the CCD layout, 
each stripe is scanned at least twice, with six interleav¬ 
ing camera columns. These two scans may have differ¬ 
ent seeing and sky noise parameters, yielding the possi¬ 
bility of a “cat scratch” pattern of imaging systemat¬ 
ics. The SDSS imaging pipeline divides each camera 
column into “fields”, which are rectangular patches of 
size 145 arcmin^. Due to the interleaving nature of the 
scanning survey, there are overlaps between neighboring 
fields. 

For SDSS DR8, to create a uniform single-pass survey 
and resolve field overlaps, each field is given a score based 
on the r-band seeing, the sky brightness in r, and an 
estim ate of photometricity of the night (Aihara et al. 
2011). Fields with higher scores are determined to be 


“primary”, with an additional percolation step that tries 
to keep regions as homogeneous as possible. However, 
as shown below there can still be a large variation in 
systematics over small scales at the edges of stripes and 
camera columns. 

At the end of this procedure, each point on the sky 
is associated with a single primary held. The primary 
window function thus determined is available from the 
SDSS-III websit^ To implement the window function, 
the sky is b roken into many mangle polygons (Swanson 
et al. 2008), each of which is linked to a primary held. 
In this way each point in the survey can be uniquely as¬ 
sociated with a point spread function (PSF) full width 
half maximum (FWHM), sky level, and sky noise in all 
hve bands. In addition, the observation time can be used 
to compute the airmass (D. Schlegel, private communi¬ 
cation) . 

For the analysis in this paper, we have con¬ 
verted the window function into HEALPix maps with 
NSIDE=2048. Each pixel has a size of 2.95arcmin^, 
which is of sufficiently high resolution to capture the hne 
scale variations of the systematics maps. However, for 
pixel-level precision then the full mangle polygon descrip¬ 
tion of the window function should be used. 


3. MEASURING SURVEY DEPTH 
3.1. Relating Survey Photometry to Survey Depth 

Let rugai be the magnitude of a galaxy. If F is the flux 
of a source in units of nanomaggies, the flux is related to 
galaxy magnitude via 

where m-zp = 22.5. If S is the number of signal photons 
reaching the detector, one has that the expectation value 
for S is related to the galaxy flux via 


^ We note that we produce our maps over the full DR8 footprint, 
but all of our verification in this paper are restricted to the BOSS 
footprint. 


(^) = kt,sF ( 2 ) 


http://data.sdss3.org/datamodel/files/PHOTO_RESOLVE 
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where fc is a proportionality constant, and teff is the effec¬ 
tive exposure time. The corresponding number of noise- 
photons N is 

(AT) = kteffFnoise (3) 

where Fnoise is the effective noise flux. The noise flux can 
have multiple origins, e.g. sky, read noise, etc. Here, we 
have lumped any such contributions into a single term. 
Given N signal photons, an unbiased estimate for the 
galaxy flux is 

F-^-F- (4) 

^ — 7 , noise* V^/ 

/Cteff 


Since both S and N follow Poisson statistics, the variance 
in F is 


2 _ {^) + {^) _ ^ Fnoise 
kHlg kteff 


( 5 ) 


The corresponding magnitude error is 
2.5 ap 
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( 6 ) 


The above expression relates the flux of a galaxy to its 
corresponding magnitude error. In practice, we recast 
^noise In terms of the lOtr limiting flux F^^ defined via 
Fiim/cTF = 10. Solving for i^noise) we have 

p . — _ p 

^ noise — 10^ ^ 

The limiting magnitude mum is simply the magnitude 
associated with a galaxy of flux Fum. 

Equations and [^together define our model for the 
magnitude error crm(m.; miini,teff) of a galaxy of magni¬ 
tude m in a survey patch of limiting magnitude mum 
and effective exposure time teff- Note that in practice, 
the effective time teff is always accompanied by the nor¬ 
malization constant fc, which can be absorbed into the 
definition of the parameter teff- 

A note on dust and dereddened magnitudes: 
Note that in the above derivation, the relevant quan¬ 
tity is the observed magnitude m, or the observed flux 
F. Whether the flux has been attenuated due to dust or 
not is completely irrelevant; the only thing that matters 
is the total arriving flux. Consequently, one should not 
use dereddened magnitudes when estimating the limiting 
magnitude of a survey. 

Conversely, we emphasize that the limiting magnitude 
we recover is a property of the survey. To relate our mea¬ 
sured limiting magnitude to a dereddened limiting mag¬ 
nitude, one must apply the necessary reddening correc¬ 
tions. As part of the data release associated with this pa¬ 
per, for convenience we include the necessary dust maps 
in the equatorial coordinate system that we use. We have 
chosen not to apply these maps to allow for the use of 
different Galactic dust models. 


3.2. Computing the Depth Maps 

As a practical application, we fit our model to SDSS 
DR8 data. Specifically, given a patch of sky, we fit our 


model of magnitude errors to the reported magnitude 
errors so as to derive the effective limiting magnitude and 
exposure times. In this procedure, we assume that the 
errors as quoted are correct. Observationally, while we 
have magnitude error estimates, we do not have “errors 
on the error”, so it is not obvious how to best fit the 
observational data to our model. In order to minimize 
the impact of gross outliers, we have chosen to fit our 
model by minimizing the total absolute deviation from 
the model, i.e. we minimize the cost function 

F(?7riim,teff) = ^ ^ l^m l^lim i^eff) | (§) 


where the sum is over all galaxies within a given sky 
patch. The funct ion is minimized usin g the downhill- 
simplex method of Nelder & Mead (1965) as implemented 
in the IDL AMOEBA function, and errors on the parameters 
are derived by bootstrap resampling the galaxies 50 times 
and refitting. 

Figurefllillustrates this technique for a single HEALPix 
pixel of NSIDE=256 in DR8. We fit all galaxies with 
signal-to-noise (S/N) greater than 5 in the pixel, which 
is deep enough to get a good measurement of the lOcr 
magnitude limit, but shallow enough to ensure that the 
SDSS “luptitudes” (arcsinh magnitudes) are equivalent 
to logarithmic magnitudes. In all panels in the figure, 
red dashed lines are our model fit. The black dotted 
lines show our recovered lOcr limiting magnitude. 

The different panels corresponds to different magni¬ 
tude definitions. The top panel uses SDSS FIBER2MAG, 
an aperture magnitude of diameter 2". We see that in 
this case the model provides an excellent description of 
the data, as we would expect. The outlier points are from 
regions that have a different local depth than the bulk 
of the relatively large pixel. Section demonstrates how 
we deal with this issue in our final depth map estimates. 
The middle and lower panels show the i-band MODEL and 
CMODEL magnitudes, which are derived from model fits to 
the galaxy data. We see that in this case there is signif¬ 
icant scatter due to the fact that the photometric noise 
for these model-fit apertures must necessarily depend on 
additional variables such as galaxy size. Nevertheless, 
it is clear that our model still provides a reasonable de¬ 
scription of the relation between galaxy magnitude and 
magnitude error, especially at the faint end where we are 
trying to estimate the limiting magnitude. 

In order to achieve a reliable fit of merr vs. mgai, we 
find we require at least 100 galaxies with signal-to-noise 
greater than 5. At the typical limiting magnitude of 
the SDSS survey, this is achieved by pixelizing the sky 


with HEALPix NSIDE=256, or 190arcmin^ per pixel. 
Note this is significantly coarser than the variations of 
the systematics maps, which introduces some noise to 
our fits which we will quantify in Section |4.4[ In some 
cases, especially at the survey edges and for most of the 
survey in the much shallower u band, there are insuffi¬ 
cient galaxies to perform the fit with NSIDE=256. When 
this happens we expand to the next coarsest pixeliza- 
tion (NSIDE=128), and if this still does not yield enough 
galaxies we use NSIDE=64. 

Figure shows the measured depth map for Tcmod- 
There is significant structure on all scales, with an am¬ 
plitude of 0.5 mag. There are several striking features 
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v^/ ^cmodelmag* * ^4 • 2 ^ 



i 

Figure 1. Magnitude error vs. magnitude for three magnitude 
definitions in a single HEALPix pixel of NSIDE=256, as follows; 
(a) SDSS FIBER2MAG, an aperture magnitude of diameter 2^^ (b) 
Same as above for MODELJIAG. (c) Same as above for CMODELJIAG. 
Red dashed lines show our model fits, and black dotted lines mark 
the recovered lOcr limiting magnitude. Note our model provides 
an excellent description of the noise for fixed metric aperture mag¬ 
nitudes, while model fitting magnitudes result in noisier fits due 
to variations introduced by galaxy size. The outlier points in the 
top panel are from regions that have a different local depth in our 
relatively large pixel. 


in the map. First, the stripe scanning pattern is imme¬ 
diately apparent, as are the “cat scratch” features from 
interleaving camera columns taken under different condi¬ 
tions. Second, when the perpendicular “iibercal” stripes 
were taken under favorable conditions, these were cho¬ 
sen as primary fields, and thus are signihcantly deeper. 
Third, the northern Galactic cap (NGP) region is typ¬ 
ically deeper than the southern Galactic cap (SGP) re¬ 
gion. The obvious exception is Stripe 82 in the south. 
This region was scanned multiple times, and thus only 
the best observations were determined to be primary 
fields for DR8. 

After performing the two-parameter fit to each of the 
pixels independently, we note that the teff and mum pa¬ 
rameters are highly correlated. This is illustrated in Fig¬ 
ure where we show a sample of pixels for both rmod 
and finod- We fit a linear model 

Inteff = a-I-6(TOiim - 21.0), (9) 

where 21.0 is a convenient pivot point. The residual scat¬ 
ter from the model in both cases is ^ 20%. 

In light of the correlation between limiting magnitude 
and survey depth, we have chosen to characterize the sur¬ 
vey depth with the single quantity miim (the more rel¬ 
evant quantity), setting the effective exposure time via 
equation Table shows the fit parameters for the 


Table 1 

teff VS. miini Fit Parameters 


Magnitude 

a 

b 

^^mod 

3.41 

1.15 

5'mod 

4.27 

0.85 

^mod 

4.53 

0.91 

^mod 

4.56 

1.00 

■^mod 

4.39 

1.34 

^cmod 

4.34 

1.70 

^cmod 

4.48 

1.23 


Note. — £eff = exp[a — 21.0)] 


magnitudes used in this paper. The efficacy of this ap¬ 
proximation is tested in the next section. 


4. RECONSTRUCTION OF THE DEPTH MAPS 
4.1. Introduction to the Problem 

In principle. Figure is exactly what we are looking 
for: a map of the limiting magnitude of the survey every¬ 
where in the sky. In practice, however, our measurement 
is severely limited by the need to have ~ 100 galaxies in a 
sky pixel in order to perform our measurement. This lim¬ 
itation renders the pixels so wide that important features 
in the small scale structure of the survey are missed. 

To overcome this difficulty, we rely on the fact that 
the observed depth map is a natural consequence of the 
survey observing conditions. For instance, better seeing 
conditions leads to deeper data. The fundamental insight 
is that the limiting magnitude at a position an¬ 

gle n is a function of systematic parameters such as sky 
noise, seeing, etc. Let then p{h) be the vector that col¬ 
lects all systematics parameters at position n. We posit 
that the function ?7iiii„(n) takes the form 

TOiim(n) = f{p{n)) (10) 


where / is some unknown function. We seek to esti¬ 
mate the function / given a series of inputs p{h) — 
the systematics evaluated for each of our coarse pixels 
— and a series of outputs ruiim — our measured depth. 
This is a problem that is well suited to machine learn¬ 
ing metho ds. Here, we re ly on the Random Forest (RF) 
technique ( Breirnjffl|2001 ), as implemented in the python 
module skleariEE 


4.2. Method 

The first step is to decide which set of systematics maps 
will be used to train the RF method. We have opted 
for an inclusive approach, and attempted to utilize all 
the systematic information that may be relevant. This 
includes the PSF FWHM, sky level, and sky sigma for 
the band being fit; t he airmass; and the value of the 
E{B — V) map from Schlegel et al. (1998). Although 
we do not expect the depth to depend on E{B — V) di¬ 
rectly, there may be some second-order correlations. As 
the model magnitude depends on the galaxy model as 
fit in the r-band, we additionally include the r-band sys¬ 
tematics maps (PSF FWHM, sky, and sky sigma) in each 
of the fits. Furthermore, the drift-scan method of SDSS 
has all the bands observed at roughly the same time, so 


® http://scikit-learn.org 
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Figure 2. Depth map for remod- The amplitude of the depth variations is 0. 5 ma g, which will be larger after correcting for Galactic 
reddening. There are numerous features that are described in the text in Section|3.2| 



the seeing and sky in the various bands are correlated 
and thus we can beat down the noise by using measure¬ 
ments from different bands. For this reason, when we 
fit the r-band depth we additionally include the i-band 
systematics maps in the fit. 

When training the RF, it is critical that the system¬ 
atics information be appropriately handled. Specifically, 
because our observed depth is measured over coarse an¬ 
gular pixels, the high resolution systematics maps must 
be averaged to the same coarser resolution employed 
in our measurements prior to training the RF. To per¬ 
form this “de-resolution” we start with the coarse target 
NSIDE=256 map. In the HEALPix nested scheme it is 
then trivial to find all the subpixels at NSIDE=2048 for 
each of the coarse pixels. We then take the average of 
the systematics value in the subpixels that are within 
the SDSS footprint, as well as tracking the fraction of 
subpixels that are outside the footprint. 

Having averaged our systematics map, we train the 
RE. In each run, we use half the pixels for training and 
half for validation. In all cases, we only use coarse pixels 
that have at least 80% of the subpixels in the footprint 
to avoid problems at the boundaries. We have also run 
checks where we train on points co-located on the sky to 
ensure that the classifier is not unintentionally picking 
up on position information encoded in nearby correlated 
pixels, and have achieved similar performance. However, 


in the case where we train on the NGP region and apply 
to the SGP region (and vice versa) the resulting perfor¬ 
mance is degraded. This hints at some subtle differences 
in the systematics properties of the north and south of 
the SDSS region. 

Finally, having trained the RF to calibrate the function 
wiim = f{p) using our low resolution systematics map, 
we apply the RF on the high resolution maps so as to 
recover a high resolution map of the limiting magnitude 

^lim(A). 

4.3. Results 

Figure illustrates the fundamental result from our 
analysis: a high resolution reconstructed map of the lim¬ 
iting magnitude For plotting purposes, the map 

was averaged to NSIDE=256. The structure of the map 
is very similar to the measured map in Figure Figure]^ 
highlights the importance of our reconstruction proce¬ 
dure. The left panel shows a zoomed-in patch of 16 deg^ 
for the low resolution depth map with ^ 200 galaxies 
per pixel. While there is a noticeable gradient across 
this patch, no higher order structure is visible. The right 
panel shows the high resolution reconstructed depth map 
for the same patch. This exhibits the characteristic “cat 
scratch” structure inherent to the interleaved scanning 
utilized in SDSS. These depth variations are tracked by 
the high resolution systematics maps, and are recovered 
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Figure 3. vs. for (a) r^od a-iid (b) *modi shown for a 

sample of pixels for clarity. Red dashed line shows a least-absolute- 
deviation fit. The residual scatter is ~ 20%. 


by our reconstruction procedure. 

Figurej^shows the residual between our observed limit¬ 
ing magnitude map, i.e., Figurej^ and our low resolution 
reconstructed map. As shown in the histogram in the in¬ 
set plot, the RMS of the residuals is Uresid = 0.03 8 ma g 
(for details on how this was computed, see Section |4.4|). 
However, there is still some structure in the residuals that 
are visible at a low level. As we will show in Section |4.4 
much of this residual is an artifact of the interplay o: 
rapidly varying systematics maps and our averaging pro¬ 
cedure. 


4.4. Estimating Noise 

We must now estimate the noise in the reconstructed 
depth maps. The most obvious estimate is the width 
of the residual distribution between the measured depth 
map and the reconstructed map as plotted in Figure 
However, this estimate suffers from some limitations. 
Most importantly, our depth map has been measured 
at a relatively coarse scale (NSIDE=256), while the aim 
of our procedure is to reconstruct the map at a fine scale 
(NSIDE=2048). As the averaging procedure adds noise, 
this will inevitably overestimate the noise in the resid¬ 
ual map as well as introduce spurious structures at field 
boundaries where the systematics vary rapidly. 

As an alternative, we use a Monte Carlo technique to 
simulate data with known systematics and repeat our en¬ 
tire analysis procedure. A comparison of the true survey 
depth map used in the Monte Carlo to the recovered sur¬ 
vey depth map will allow us to estimate the noise in the 
reconstructed maps. For our Monte Carlo simulation, we 
define the survey depth to be our estimated high resolu¬ 
tion survey depth map, and for our input galaxy catalog 
we use the SDSS galaxy catalog, defining the observed 


magnitudes as truth for the Monte Carlo simulation. Our 
procedure is as follows: 

1. Take the existing galaxy catalog, and set each (ob¬ 
served) magnitude as the “true” magnitude in each 
band. 

2. Take the high resolution reconstructed depth map 
for each band as the “true” limiting magnitude 
(and associated teff) at every point in the survey. 

3. Given the true magnitude, miim, and teS, compute 
the associated magnitude error from Eqn. Per¬ 
turb the true magnitude according to this error. 

4. Compute the dep th m ap with NSIDE=256 as de¬ 
scribed in Section [321 

5. Refit the depth map and reconstruct a simulated 
high resolution depth map as described in Sec¬ 
tion IH 

In the end, we have a high resolution reconstructed depth 
map which can be directly compared to the “true” input 
map. 

We find that the residual map between the raw mea¬ 
sured depth map and our recovered high resolution depth 
map (degraded to the resolution of the measured map) 
has a striking similarity to that of Figure Since we 
now have the high resolution “truth” in hand, we can 
determine the origin of the structures shown in Fig¬ 
ure!^ Figure shows a particularly problematic patch 
at RA = 225°, DEC = 42°. The left panel shows the 
NSIDE=256 residuals between the measured depth map 
and the averaged reconstructed map for the simulated 
data. The red streak shows a large offset of > 0.1 mag, 
very similar to that observed in Figure for the real 
data. The right panel shows the NSIDE=2048 high res¬ 
olution residuals between the reconstructed simulated 
depth map and the “true” map used as an input. We see 
that the pattern essentially disappears. Evidently, the re¬ 
gions with large coherent residuals observed in Eigure|^ 
are artifacts of the averaging process. Most importantly, 
our high resolution maps are successfully recovering high 
resolution features that are completely unresolved in the 
observed depth maps. 

In Table we summarize the noise estimates for the 
reconstructed depth maps for two different methods. The 
first, CTresidj measures the width of the residual distribu¬ 
tion between the measured depth map and the averaged 
reconstructed depth map (NSIDE=256), as in the in¬ 
set plot of Figure!^ The second, (Tsim, is the width of 
the residual distribution between the simulated recon¬ 
structed depth map and the “true” depth map. In all 
cases the width is estimated using the robust Gaussian 
histogram-fitting code histogauss.pro from the IDL 
Astronomy Librarjj^ Due to the fact that our simulated 
data cannot take into account all the various sources of 
unmodeled systematics, we consider (Jsim and CTresid to be 
lower and upper bounds on the true noise of the recon¬ 
structed depth map. 

http://idlastro.gsfc.nasa.gov/ 
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Figure 4. Reconstructed map of rcmod depth usin g RF . The map was reconstructed at high resolution (HEALPix NSIDE=2048) and 
then averaged to NSIDE=256 as described in Section [4^ for plotting and comparison to the measured map. 
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Figure 5. Zoom in of 16deg^ of depth map. Left: Low resolution (NSIDE=256) with ~ 200 galaxies per pixel used for the depth 
estimation. There is a noticeable gradient, but no higher order structure is visible. Right: High resolution (NSIDE=2048) reconstructed 
depth map showing the “cat scratches” characteristic of the interleaved scanning. These depth variations are tracked by the high resolution 
systematics maps, and are recovered by our reconstruction procedure. 
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Figure 6. Residuals of reconstructed rcmod depth compared to the observed limiting magnitude map. The histogram in the inset plot 
shows the RMS of the residuals, with a G aussian fit (red) with o-j-esid = 0.038 mag. There is some structure in the residuals that are visible 
at a low level. As shown in Section |4.4| the large outliers are caused by the averaging process used to compare the maps, but are not 
present in the underlying high resolution map. 



Figure 7. Zoom in of simulated depth map residuals. Left: Low resolution (NSIDE=256) residuals between measured d^th map 
and averaged reconstructed depth map for the simulated data. The red streak shows a large offset, similar to that of Figure Im Right: 
High resolution (NSIDE=2048) residuals. The large offset is no longer present, highlighting that our procedure is recovering the limiting 
magnitude at fine scales. 
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Table 2 

Observed and Simulated 
Scatter 


Magnitude 

^resid 

^sim 

'^mod 

0.128 

0.068 

9mod 

0.038 

0.020 

^mod 

0.036 

0.018 

^mod 

0.040 

0.020 

^mod 

0.063 

0.034 

^cmod 

0.038 

0.018 

^cmod 

0.051 

0.023 


4.5. Using 5cr Point-Source Depth 

Many sky surveys state depth in terms of 5cr point- 
source depth, and this is typically used for telescope 
exposure-time calculators. This has the advantage of be¬ 
ing relatively straightforward to compute, but as we show 
in this section this does not scale trivially to the galaxy 
limiting magnitude which we are interested in comput¬ 
ing. 

The BOSS survey have independently modeled the 5cr 
point-source depth based on the formal errors from PSF 
photometry on stellar sources (D. Schlegel, private com¬ 
munication). BOSS uses a simple model of seeing, sky 
brightness, and airmass: 

F(5a) = (11) 

where fi is the FWHM, Si is the sky flux, and 
rui is the airmass. The normalization term a = 
{0.387,0.218,0.241,0.297,0.665} and the airmass term 
k = {0.49,0.17,0.10,0.06,0.06} for bands i = {u,g,r,i,z} 
respectively. This formula predicts the 5a point-source 
depth to an accuracy of 12%, 4%, 4%, 5%, and 7%, nom¬ 
inally comparable to the accuracy of the galaxy depth 
reconstruction in this work, but for a different type of 
object. 

To compare to the galaxy depth, we calculate the 5a 
point-source depth at high resolution and then average 
the map to NSIDE=256 as above. We then apply a con¬ 
stant offset to account for the fact that the lOcr galaxy 
depth is considerably shallower than that for these point 
sources. In the case of Tcmod, this offset is 1.43 mag. 
Figure shows the residual between our observed lim¬ 
iting magnitude map (Figure]^ and the low resolution 
map reconstructed from the point source depth. The 
RMS scatter is 11%, a factor of three larger than seen 
for our galaxy depth reconstruction technique in this pa¬ 
per. More importantly, there is significant structure in 
the residuals. In Section]^ we show the impact on corre¬ 
lation function measurements of using this map instead 
of a galaxy-appropriate map. 


5. GALAXY CATALOG COMPLETENESS 

With our depth maps in hand, we can now properly 
estimate the completeness of the galaxy detection as a 
function of local depth. We make this estimation em¬ 
pirically, using two overlapping surveys that are signifi¬ 
cantly deeper than single-pass SDSS DR 8 imaging. The 


first is the SPSS Stripe 82 (“S82”) coadd catalog (An- 


nis et al. 
our DR 8 


2011), which covers ~ 205 deg overlapping 
ootprint, with limiting magnitudes ~ 2 mag 


deeper than the sin gle-pass imaging. The second is the 


CFHTLens catalog (Heymans et al. 2012 Erben et al. 


2013), whose W1 and W4 fields cover 
limiting magnitudes 


75 deg , with 
3 mag deeper than DR 8 . Along 


with being based on two different surveys, these catalogs 
sample different parts of sky and different systematics, 
and as such they are an adequate independent test of 
the detection completeness. 


5.1. Measuring Completeness 

We investigate how survey completeness depends on 
survey depth. Our method for measuring the complete¬ 
ness of galaxy detection is quite straightforward. We 
start with the “truth,” which is the designation we give to 
the deeper reference catalog (either S82 or CEHTLens). 
Our procedure is then: 

1. Apply the BOSS mask to both surveys. 

2. Pixelize the footprint with HEALPix using 
NSIDE=2048. 

3. Bin the pixels according to local depth, with 10 bins 
over S82 and 5 bins over GFHTLens (because of 
the smaller area). Note that these pixels need not 
be contiguous; rather, they must share a common 
local depth. 

4. For each depth bin, compute a median zero-point 
correction to put galaxies on a “DR 8 ” scale. This 
removes any overall shifts due to filter and magni¬ 
tude measurement method. 


5. Measure the detection rate of galaxies as a function 
of true magnitude in each depth bin. This is the 
desired completeness. 

6 . Bootstrap resample the pixels in each depth bin 
1000 times, remeasuring the depth each time in 
order to determine observational uncertainties. 


7. Fit for the completeness as a function of magni¬ 
tude using the functional form detailed below. The 
bootstrap resampling of the input pixels is used to 
estimate errors on parameters. 


Figure shows the completeness estimated for rcmod 
detections for the median DR 8 depth in the S82 fields. 
We model the completeness with a simple error function 
model: 


c = (e/ 2 ) 


1 — erf 


m - TO 50 


( 12 ) 


where e is the overall efficiency of detection at the bright 
end, TO 50 is the magnitude at which the completeness is 
50%, and w is the (Gaussian) width of the rollover. The 
red dashed line in the plot shows the best-fit model, and 
the blue dotted line shows the median lOcr limit in the 
sampled region. The error bars are smaller than the data 
points at all but the brightest magnitudes. Although the 
model fit is not perfect, it clearly provides a reasonable 
description of the shape of the completeness function. 
The end result is that we now have measurements of mso 
and w in bins of different survey depth rniim- 
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Figure 8. Residuals of reconstructed rj-mod depth using the formula for the 5cr point-source depth, compared to the measured depth map. 
The RMS scatter is 11%, a factor of three larger than seen for our galaxy depth reconstruction technique in this paper. More importantly, 
there is significant structure in the residuals. 




Figure 9. Completeness estimated for rcmod for Stripe 82 data 
for the median DR8 depth in the Stripe 82 fields. R ed d ashed line 
shows the fit of the functional form given by Eqn. |12| The blue 
dotted line shows the lOcr limiting magnitude. The error bars are 
smaller than the data points at all but the brightest magnitudes. 


5.2. Completeness Results 

With these measurements in hand we determine how 
the galaxy completeness parameter m^o depends on the 
survey depth parameter mum. Figures 10 and 11 show 


the completeness as a function of the local limiting mag¬ 
nitude measured for rcmod and tcmod respectively. Blue 
circles show the results for S82 regions, and red squares 
for CFHTLens. The S82 region s te nd to be deeper than 
average (as discussed in Section 3.2), and thus we sample 
a wide range of local depth with some overlap between 
the two data sets. The two data sets give consistent re¬ 
sults, with a small systematic offset of ~ 0.01 magnitude. 
Solid error bars are those estimated from the bootstrap 
resampling, and the dashed error bars are an estimate of 
the systematic error required to yield x^/dof = 1 for the 
linear ht. 

Table 1^ shows the final completeness parameter hts for 
TO 50 anaw. We note that both 77159 and w are strong 
functions of the local limiting magnitude. Thus, in order 
to estimate the completeness at any point in the survey, 
we must first know the local depth. 


6 . IMPACT ON CORRELATION FUNCTIONS 

As an example of the use of these depth maps, we 
now investigate the impact on the estimation of correla¬ 
tion functions. The measurement of the correlation func¬ 
tion w{9) requires a set of random poin ts that describes 
the w indow function of the observable (Tandy & Szalay 
1993|). To first order, these random points can be sam¬ 
pled uniformly from the geometric mask of the survey. 
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Table 3 

Completeness Fit Parameters 


Magnitude 

Q-comp 

^comp 

a^j 


^cmod 

^cmod 

21.67 ± 0.02 
21.27 ± 0.01 

0.53 ±0.04 
0.41 ± 0.02 

0.24 ± 0.01 
0.206 ±0.003 

-0.26 ±0.03 
-0.16 ±0.02 


Note. — mso = flcomp + ^>comp(lim — 21) ] w = a^n + foi(;(lim - 21) 


22.00 

21.95 

21.90 
I' 21.85 
21.80 

21.75 


Figure 10. Main panel: Completeness parameter mso, the mag¬ 
nitude at which 50% of the galaxies are detected, as a function 
of local depth measured from rcmod- Blue circles are from S82, 
red squares from CFHTLens. Black dashed line is a linear fit. In- 
se£:Width parameter w as a function of local depth. Symbols are 
the same as the main panel. 
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- ■ CFHTLens ^ - 
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Figure 11. Same as Figure [To| with limiting magnitude measured 
from icmod- 


Although this works well at brighter magnitudes, when 
pushing the limits of a photometric survey this is certain 
to break down. In this section, we measure a simulated 
correlation function while comparing different methods 
of generating random points to show the utility of our 
method at faint magnitudes. 

6.1. Generating a Simulated Survey 

Our first task is to generate a simulated DR8-like sur¬ 
vey with known systematics and known correlation func¬ 
tion. We follow the same general plan as in Section [4^ 
using a Monte Carlo technique to simulate data with uni¬ 
form density (zero correlation). Consequently, any cor¬ 


relations observed in the data must necessarily be due to 
systematics being imprinted into the galaxy catalog. 

There are two primary differences between this simu¬ 
lation and that from Section lT4l The first is that we are 
building in the completeness estimation from Section 
The second, is that we sample galaxies by color as wml 
as magnitude to enable color selection in random point 
generation. Our procedure is as follows: 

1. Generate simulated galaxies uniformly with zero 
correlation and approximately the same density as 
DR8 with icmod < 21. 

2. Give each galaxy a true magnitude vector sampled 
from the (deeper) S82 catalog. 

3. Estimate the completeness for each galaxy based 
on icmod- Remove galaxies randomly according to 
the completeness estimate. 

4. For each magnitude of each of the remaining galax¬ 
ies, perturb the input magnitudes in accordance to 
the input survey de pth map to yield a simulated 
survey as in Section [T^ 

As before, given this simulated survey we recompute 
the depth maps, fit the depth map, reconstruct a sim¬ 
ulated high resolution depth map, and estimate the 
completeness from the simulated data. It is these re- 
estimated maps that are used in the following tests. 


6.2. Generating Random Points 

In order to measure the correlatio n function in real- 
space, we make use of the C0RR2 code (Jarvis et al.|2004| 
for computing the Landy-Szalay statistic: 


n = 


DD - 2DR -b RR 
RR ■ 


(13) 


For our measurements, D refers to the simulated galax¬ 
ies, and R refers to the random points used to sample the 
survey mask. We have chosen to compare the impact 
of four different methods of generating random points. 
These are: 


1. “Uniform”: Only the geometric mask has been ap¬ 
plied. 

2. “Low Resolution”: The low resolution measured 
depth maps (NSIDE=256) are used. We also apply 
the completeness correction. 

3. “Reconstructed”: The high resolution recon¬ 
structed depth maps (NSIDE=2048) are used. We 
also apply the completeness correction. 
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4. “PSFMAG”: The 5cr point-source psf magnitude 
depths from Section 4^ are used, along with the 
derived constant offsets so the values match the 
lOcr galaxy depths on average. We also apply the 
completeness correction. 


In all cases, the generation of random points follows 
the method for building the simulated survey galaxies in 
the first place. The primary difference is that we use the 
new reconstructed depth maps and re-estimated com¬ 
pleteness. Therefore, the difference between the simu¬ 
lated survey galaxies and the random points is the noise 
in the estimation of the depth and completeness maps. 


6.3. Results 

In the interest of simplicity, we measure the correla¬ 
tion function for two magnitude cuts in i-band, and ig¬ 
nore the impact of the color selection. The first mag¬ 
nitude cut estimates the impact of depth variations on 
the correlation function for the faintest BOSS CMASS 


galaxies (19.8 < icmnd < 2 0.0; White et al. 2011 Ross 
et al. 2011 Ho et al. 2012). The second magnitude cut 
estimates the impact of depth variations on photometri¬ 
cally detected galaxies near the average lOcr survey limit, 
20.9 < icmod < 21.0. 

Figure [T^ shows the absolute value of the two-point 
correlation function |rc(0)| estimated for various meth¬ 
ods of generating random points. The dotted bl a ck lin e 
shows the typical error on w{d) from |Ho et al. (2012). 
Any additive bias on the correlation function below this 
level is sub-dominant. In the top panel, for the faintest 
CMASS galaxies with 19.8 < icmod < 20.0, we can see 
that uniform randoms (sold black line) are adequate to 
compute w{0). Similarly, the other random point genera¬ 
tion schemes, with the same geometric mask, all perform 
adequately. The only possible exception is the “Low Res¬ 
olution” method (magenta dotted line) which may im¬ 
part some structure at very small scales. In the bottom 
panel, for the galaxies near the lOcr photometric limit 
of the survey (20.9 < fcmod < 21.0), the choice of ran¬ 
dom points is much more significant. Using uniform ran¬ 
dom points (with the geometric mask) leads to significant 
residuals at all scales. The “PSFMAG” estimate of the 
depth (blue dot-dashed line) is better than uniform, but 
still has problems at all scales. The “Low Resolution” 
method works at large scales but not at small scales, 
while our full “Reconstructed” method (red dashed line) 
yields the best results, with only a small bias at the small¬ 
est scales. Therefore, we are confident that this method 
accurately captures depth and completeness variations 
down to faint magnitudes. 


7. RESULTS AND SUMMARY 

In this paper, we have introduced a new way of empir¬ 
ically constructing depth maps from large photometric 
survey data. By combining low resolution depth maps 
with high resolution maps of survey systematics, we can 
use machine learning to reconstruct high resolution depth 
maps for typical galaxies measured using any arbitrary 
magnitude definition. Using SDSS DR8 as an example, 
our method can reconstruct the depth at an arbitrary 
point in the survey to within ~ 2%, thus making it possi¬ 
ble to accurately measure correlation functions for galax¬ 
ies detected near the limiting magnitude of the survey. 



Figure 12. Absolute value of the correlation function |ui(6)| for 
our simulated survey data using different sets of random points. 
Top Panel: Magnitude selection equivalent to the faintest BOSS 
CMASS galaxies, 19.8 < < 20.0. The bl ack dotted line shows 

the typical error on w(9) from Ho et al.|||2012ll. The solid black line 
shows the bias from using unitorm randoms (with the geometric 
mask). The magenta dotted line shows the bias from the “Low 
Resolution” randoms. The red dashed line shows the bias from the 
full “Reconstructed” randoms, and the blue dot-dashed line shows 
the bias from the “PSFMAG” randoms. Bottom Panel: Same as 
above, for galaxies near the lOtr survey limit, 20.9 < fcmod < 21.0. 
The impact on w{9) for uniform randoms is large at all scales. Only 
the full reconstructed randoms show nearly unbiased results at all 
scales. 


We emphasize that this method requires the knowledge 
of both a precise geometric mask, as well as accurate 
maps of the various systematics that affect galaxy pho¬ 
tometry, especially seeing and sky brightness. Further¬ 
more, the reported photometric errors in the data must 
be correct: any biases in the estimated photometric un¬ 
certainties of the detected galaxies will propagate into 
the survey depth maps. 

In order to accurately build sets of random points suit¬ 
able for correlation function measurements, a depth map 
is insufficient. We must also estimate the completeness 
function for galaxy detection. To estimate this func¬ 
tion directly from the data, it is necessary to have data 
that is significantly deeper than the main survey. Fortu¬ 
nately, this is not only available for SDSS in the Stripe 82 
coadds, but it is a common aspect of current and future 
wide-held surveys, inclu ding the Medium- Deep Survey 
helds in Pan-STARRsl (Kaiser et al. 2002) and the su¬ 


perno va search helds in DES (I’he DES Collaboration 


20051. We have shown that the completeness can be 


parametrized by a simple function for SDSS, and that 
the parameters of this function depend on the local sur¬ 
vey depth. In this way, the depth map reconstruction 
is an essential prerequisite for estimating the local com¬ 
pleteness of the survey. 

As a sample of the utility of these depth maps, we 
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demonstrate a new way of building random points for 
correlation function measurements that incorporates the 
full knowledge of the survey depth and local complete¬ 
ness. Uniform random points are sufficient when one is 
far from the photometric limit of the survey, as is the 
case for BOSS CMASS galaxies. These galaxies are rela¬ 
tively unaffected by variations in the local depth for two 
main reasons. First, the drift-scan technique ensures that 
SDSS is a relatively uniform survey, without large varia¬ 
tions in depth. Second, we note that while SDSS uses the 
same telescope for both photometric and spectroscopic 
measurements, by necessity any spectroscopic target se¬ 
lection (such as CMASS) will necessarily be significantly 
brighter than the photometric limit for any reasonable 
balance of exposure time between the two surveys. Thus, 
it is not surprising that for the CMASS sample, designed 
for spectroscopic follow-up, uniform random points (with 
a proper geometric mask) are sufficient. 

However, for ongoing and future photometric sur veys, 


this will not be the case. In particul ar, both DES ( The 
DES Collaboration 20051 and LSST (Ivezic et al. 2008| ) 


are purely photometric surveys where the best science 
will not be achieved by limiting measurements to galax¬ 
ies more than a magnitude brighter than the survey 
limit. Furthermore, with tiled observations and chip 
gaps, depth variations will be much more significant than 
in SDSS, especially in the early phases before complete 
coverage is achieved. In addition, surveys such as the 
DECam Legacy Survey (DECaL^^ will be used for tar¬ 
get selection for the Dark Energy spectrographic Instru¬ 
ment Survey (DESf^ are not significantly deeper than 
the spectroscopic targeting. Therefore, we expect this 
method of building depth maps will be useful beyond 
the current application to SDSS. 
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APPENDIX 

A. RECONSTRUCTED DEPTH MAPS 

The maps are available at http://risa.stanford.edu/redmapper, The list of maps are in Table[^ Each of them 
is in HEALPix FITS format, NSIDE=2048, ring ordered, equatorial coordinates. 


http://legacysurvey.org 
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14 


Rykoff et al. 


Table 4 

List of HEALPix Maps 


Band Type 


Name 


u 

9 

r 

i 

z 

r 

i 


model 

model 

model 

model 

model 

cmodel 

cmodel 

E{B-V) 


sdss_dr8_nodered_nside2048_u_model_10sigma. f its .gz 
sdss_dr8_nodered_nside2048_g_model_10sigma. f its .gz 
sdss_dr8_nodered_nside2048_r_model_10sigma. f its .gz 
sdss_dr8_nodered_nside2048_i_model_10sigma. f its .gz 
sdss_dr8_nodered_nside2048_z_model_10sigma. f its .gz 
sdss_dr8_nodered_nside2048_r_cinodel_10sigina. f its .gz 
sdss_dr8_nodered_nside2048_i_cinodel_10sigina. f its .gz 
sdss_dr8_ebv_sfd98jnside2048 .f its .gz 


Note. — All files are in HEALPix PITS format, NSIDE=2048, ring ordered, 
equatorial coordinates. 



